
options(rasterTmpDir='C:/temp/')
cott_fn <- 'V:/00_GLB/001_farming/SPAM_2010_v2r0/spam2010V2r0_global_P_COTT_A.tif'

##
## IMPORTANT PATH to Rtools to run Cpp
##
Sys.setenv(PATH = paste("C://WBG//Utility//Rtools//bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "C://WBG//Utility//Rtools//mingw_$(WIN)//bin//")

require(Rcpp)

output_version <- '2020-09-22' # Sys.Date()

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones)

##
##
## Cotton : SPAM 2010v2
##

output_version <- '2020-10-20'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)


cott.glb <- raster(cott_fn)
# manually set outside for 5x5
ntl.extent <- c(-5,30,-5,30)
cott <- crop(cott.glb,ntl.extent,snap='out')

## - create measure avg same cell, avg 5 cells from cell, avg 10 cells from cell
cott.diam.crop <-  sqrt(raster::area(cott))
dd050km = 50 / cellStats(cott.diam.crop,'mean')
lendd050km = res(cott)[1] * dd050km

dd100km = 100 / cellStats(cott.diam.crop,'mean')
# 1deg = 111 km  (or 60 nautical miles)
lendd100km = res(cott.diam.crop)[1] * decdeg100km

fw5x5 <- focalWeight(cott, lendd050km, type='circle')
fw10x10 <- focalWeight(cott, lendd100km, type='circle')

cott050 <- focal(cott, w=fw5x5, fun=sum,na.rm=TRUE)
cott100 <- focal(cott, w=fw10x10, fun=sum,na.rm=TRUE)

ntl.zones$cott <- exactextractr::exact_extract(cott,ntl.zones,fun='sum')
ntl.zones$cott050 <- exactextractr::exact_extract(cott050,ntl.zones,fun='sum')
ntl.zones$cott100 <- exactextractr::exact_extract(cott100,ntl.zones,fun='sum')

output_version <- '2020-10-20'

save.dta13(st_drop_geometry(subset(ntl.zones,select=c(OBJECTID,cott,cott050,cott100))), 
           paste0(wkdir,"/proc_data/COTT",output_version,".dta"))
